# polygonclass by Kristoffer Josefsson

import math
import string
import random

from OpenGL.GL import *
#from scipy import *
#from Numeric import *

EPSILON = 10e-6

class vector(list):
	#def __init__(self,v):
	#	self.v = v
	#def __repr__(self):
	#	return `self`
	def __add__(self,v):
		return vector(map(lambda x,y: x+y, self, v))
	def __sub__(self,v):
		return vector(map(lambda x,y: x-y, self, v))
	
	def __rmul__(self,s):
		self = vector(map(lambda x: x*s, self))
		return self
	def __imul__(self,s):
		self = vector(map(lambda x: x*s, self))
		return self
	def __div__(self, s):
		sinv = 1/s
		self = vector(map(lambda x: x*sinv, self))
		return self

	def atov(self,a):
		return vector(map(lambda x: vector(x), a))
	
	def dot(self,v):
		return sum(map(lambda x,y: x*y, self, v))
	# dot product as * operator
	def __mul__(self,v):
		return self.dot(v)

	def __abs__(self):
		return math.sqrt(self*self)

	def cross(self,v):
		if (self.dim() == v.dim()) and (self.dim() == 3):
			a = self
			b = v
			return vector([a[1]*b[2]-a[2]*b[1],a[2]*b[0]-a[0]*b[2],a[0]*b[1]-a[1]*b[0]])
		else:
			raise NameError, 'vector.cross only supports 3-vectors'

	# crossproduct as % operator
	def __mod__(self,v):
		return self.cross(v)
	
	def normalize(self):
		len = abs(self)
		if len != 0:
			self=[x/len for x in self]
		return vector(self)

	def dim(self):
		return len(self)

class polygon:
	
	def __init__(self,filename=""):
		self.vert = [];
		if filename == "":
		#	self.compute()
			self.makeTorusKnot()
		else:
			self.read(filename)

		self.holComp = True

	def __str__(self):
		return('Vertices:\t'+`len(self.vert)`+'\n'
				'Dimension:\t'+`self.vert[0].dim()`+'\n'
				'Closed:\t'+`self.closed`+'\n'
				'Total Curvature:\t'+`self.tc()`+'\n'
				'Holonomy:\t'+`self.holonomy()`)+'\n'

	def update(self, a):
		return

	def tc(self):
		return sum(map(abs,self.curv))

	def holonomy(self):
	#	return (sum(self.tors))%(2 * math.pi)
		return abs(math.acos(self.norm[0]*self.norm[-1]))

	def read(self,filename):
		self.vert = []
		if(filename != ""):
			f = open(filename)
			for s in f.readlines():
				tlist = string.split(s,"\t")
				new = vector(map(float,tlist))
				self.vert.append(new)

			self.closed = False;
			if self.vert[0] == self.vert[-1]:
				# closed curve
				self.closed = True
				self.vert.pop()
			else:
				self.closed = False
				
			self.prepare()

	# these two from peters example
	def trail(self,func, s0=0, s1=2*math.pi, N=64, closed=1):
		return [vector(func(s0+(s1-s0)*i/N)) for i in range(N)]
	
	def torusKnot(self,r=1., n=2, k=3, N=64):
		return self.trail(lambda s: [r*math.cos(n*s)*(2+math.cos(k*s)),
								r*math.sin(k*s),
								r*math.sin(n*s)*(2+math.cos(k*s))], N=N)

	def makeTorusKnot(self):
		r = random.randint(1,10)
		n = random.randint(1,10)
		k = random.randint(1,10)
		N = 2**random.randint(5,8)
		self.closed = True
		self.vert = self.torusKnot(r,n,k,N)
		self.prepare()

	def compute(self):
		if self.vert == []:
			self.makeTorusKnot()
			
	def randomize(self):
		return

	# calculate the Bishop frame for the curve etc.
	def prepare(self):

		# we need tangents
		self.tan = []
		for i in range(len(self.vert)):
			self.tan.append((self.vert[(i+1) % len(self.vert)] - self.vert[i]).normalize())

		# curvature at each vertex
		# question, how get rid of loop
		self.curv = []
		for i in range(len(self.tan)):
			v1 = self.tan[i]
			v2 = self.tan[(i+1) % len(self.tan)]
			ainv = v1*v2
			self.curv.append((math.acos(ainv)))

		# bishop normals
		if self.dim() >= 3:
			self.norm = []
			# first normal (arbitrary orthogonal to tangent)
			# pick vector in plane to build plane containing tangent
			v = vector([1,2,3]) # small enough odds this is paralell ?
			n = (v%self.tan[0]).normalize()

			self.norm.append(n)
			
			# now project this normal onto the next normalplane
			# and find the bishop normal
			for i in range(1,len(self.tan)+1):
				#kinv = 1/self.curv[i]
				#n = (self.tan[(i+1) % len(self.tan)] - self.tan[(i-0) % len(self.tan)]).normalize()
				#self.norm.append(n)

				# plane equation is self.tan(i+1)
				# direction to project from is self.tan(i+1)-self.tan(i)
				#m = self.tan[(i+1) % len(self.tan)] - self.tan[i]
				#n = m%(self.tan[(i+1) % len(self.tan)]%self.tan[i])
				#self.norm.append(n.normalize())

				# now then...
				v = self.tan[i-1]
				n0 = self.norm[i-1]	# gammal normal
				N = (v%self.tan[i%len(self.tan)])%(v-self.tan[i%len(self.tan)])	#plan
				a = v.dot(N)
				if abs(a) > EPSILON:
					p = n0 - ((n0.dot(N))/a)*v
					p = p - (p.dot(v))*N
				else:
					p = vector((1,0,0))
				self.norm.append(p.normalize())

			self.norm.pop(0)

			# binormals
			self.binorm = []
			for i in range(len(self.tan)):
				b = self.tan[i].cross(self.norm[i])
				self.binorm.append(b.normalize())

			# normals again
			#self.norm = []
			#for i in range(len(self.binorm)):
			#	n = self.tan[i].cross(self.binorm[i])
			#	self.norm.append(n.normalize())

			# torsion
			#self.tors = []
			#for i in range(len(self.norm)):
			#	tinv = (self.norm[(i+1) % len(self.norm)]*self.norm[i]) - 1e-10
			#	ang = math.acos(tinv)
			#	self.tors.append(ang)

	def dim(self):
		return self.vert[0].dim()

	def append(self,v):
		self.vert.append(v)
		self.prepare()

	def subdivide(self,n, a = 9./8.):
		npts = len(self.vert)
		dim = self.vert.dim()

		sdmatrix = array()

		array([[0,1,0,0,0], [(1-a)/2,a/2,a/2,(1-a)/2,0], [0,0,1,0,0],[0,(1-a)/2,a/ 2,a/2,(1-a)/2],[0,0,0,1,0]], Float)


class polyGL(polygon):

	def __init__(self,filename=""):
		polygon.__init__(self,filename)
		self.listids = glGenLists(10)
		self.redraw = True

	def draw(self):


		glPushAttrib(GL_ENABLE_BIT | GL_LIGHTING_BIT)
		glDisable(GL_LIGHTING)
		glColor3f(1.0,0.0,0.0)
		glLineWidth(2)
		glBegin(GL_LINE_STRIP)
		for v in self.vert:
			glVertex3fv(v)

		if self.closed:
			v = self.vert[0]
			glVertex3fv(v)
		glEnd()
		glPopAttrib()


	def drawQuad(self):
		glBegin(GL_QUADS)
		for v in self.vert:
			glVertex3fv(v)

	#	if self.closed:
		if 0:
			v = self.vert[0]
			glVertex3fv(v)
		glEnd()
		
	

##	old and unupdated so won't use them
##	def drawTan(self):
##		glColor3f(1.0,0,0)
##		glLineWidth(2)
##		glBegin(GL_LINES)
##		for i in range(len(self.tan)):
##			v = self.vert[i]
##			t = self.tan[i]
##			vt = v+t
##			glVertex3fv(v)
##			glVertex3fv(vt)
##		glEnd()
##
##	def drawNorm(self):
##
##		glColor3f(0,1.0,0)
##		glLineWidth(2)
##		glBegin(GL_LINES)
##		for i in range(len(self.vert)):
##			v = self.vert[i]
##			t = self.norm[i]
##			vt = v+t
##			glVertex3fv(v)
##			glVertex3fv(vt)
##		glEnd()
##
##	def drawBinorm(self):
##
##		glLineWidth(2)
##		glBegin(GL_LINES)
##		for i in range(len(self.binorm)):
##			glColor3f(0,0,self.curv[i])
##			v = self.vert[i]
##			t = self.binorm[i]
##			vt = v+t
##			glVertex3fv(v)
##			glVertex3fv(vt)
##		glEnd()

	def drawHol(self):

		glPushAttrib(GL_ENABLE_BIT | GL_LIGHTING_BIT )
		glDisable(GL_LIGHTING)
		#glDisable(GL_DEPTH_TEST)
		glLineWidth(10)
		glBegin(GL_LINES)

		glColor4f(0,1.0,0,1.0)
		v = self.vert[0]
		t = self.norm[0]
		vt = v+t
		glVertex3fv(v)
		glVertex3fv(vt)

		if not self.holComp:
			if self.holonomy() > 1e-5 or self.holonomy() < 1e-5:
				glColor3f(1.0,0,0)
				v = self.vert[0]
				t = self.norm[-1]
				vt = v+t
				glVertex3fv(v)
				glVertex3fv(vt)

		glEnd()
		glPopAttrib()


	def makeFrameCircle(self,norm,binorm,r,s,o=0):
		c = []
		for i in range(s):
			c.append(r*math.cos(o + (i * 2 * math.pi / s))*norm + r*math.sin(o + (i * 2 * math.pi / s))*binorm)
		return c

	def drawTube(self, radius = 0.5, segments = 16):


		#glPolygonMode(GL_FRONT_AND_BACK,GL_LINE)

		glPushAttrib(GL_ENABLE_BIT | GL_LIGHTING_BIT)
		glEnable(GL_LIGHTING)

		glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, [0.7, 0.6, 0.7, 0.1])
		glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, [0.9, 0.9, 0.7, 0.1])
		glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 20.0)
		glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT, [0.5, 0.7, 0.1, 0.1])
		glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, [0.6, 0.8, 0.5, 0.1])

		glBegin(GL_QUAD_STRIP)

		c0 = c2 = self.makeFrameCircle(self.norm[0],self.binorm[0],radius,segments)
		for i in range(len(self.vert)-1):
			# this can be better, for example consider length of segment
			# then must be moved further into loop
			if self.holComp == True:
				offset = i*self.holonomy()/len(self.vert)
			else:
				offset = 0

			c0p = c1 = c2
			c2 = self.makeFrameCircle(self.norm[i+1],self.binorm[i+1],radius,segments,offset)
			for s in range(segments):
				v1 = self.vert[i]+c1[s]
				v2 = self.vert[i+1]+c2[s]
				glNormal3fv(c1[s])
				glVertex3fv(v1)
				glVertex3fv(v2)

			v1 = self.vert[i]+c0p[0]
			v2 = self.vert[i+1]+c1[0]
			glNormal3fv(c0p[0])
			glVertex3fv(v1)
			glVertex3fv(v2)

		if self.closed:
			for s in range(segments):
				v1 = self.vert[-1]+c2[s]
				v2 = self.vert[0]+c0[s]
				glNormal3fv(c2[s])
				glVertex3fv(v1)
				glVertex3fv(v2)

			v1 = self.vert[-1]+c0p[0]
			v2 = self.vert[0]+c1[0]
			glNormal3fv(c0p[0])
			glVertex3fv(v1)
			glVertex3fv(v2)

		glEnd()

		glPopAttrib()


	def drawBand(self,width=0.5):


		glPushAttrib(GL_ENABLE_BIT | GL_LIGHTING_BIT)
		glEnable(GL_LIGHTING)

		glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, [0.4, 0.6, 0.7, 0.5])
		glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, [0.2, 0.9, 0.7, 0.5])
		glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 10.0)
		glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT, [0.3, 0.2, 0.1, 0.5])
		glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, [0.2, 0.1, 0.1, 0.7])

		#glPolygonMode(GL_FRONT_AND_BACK,GL_LINE)

		w = width-1e-8

		glBegin(GL_QUAD_STRIP)

		for i in range(len(self.vert)):

			if self.holComp == True:
				offset = i*self.holonomy()/len(self.vert)
				n1 = math.cos(offset)*self.norm[i] + math.sin(offset)*self.binorm[i]
			else:
				n1 = self.norm[i]

			v1 = self.vert[i% len(self.vert)]
			v2 = self.vert[i% len(self.vert)]+w*n1

			glNormal3fv(n1)
			glVertex3fv(v1)
			glVertex3fv(v2)

		if self.holComp == True:
			offset = self.holonomy()
			n1 = math.cos(offset)*self.norm[-1] + math.sin(offset)*self.binorm[-1]
		else:
			n1 = self.norm[-1]

		if self.closed:
			glNormal3fv(n1)
			glVertex3fv(self.vert[0])
			glVertex3fv(self.vert[0]+w*n1)

		glEnd()

		glPopAttrib()


	def drawCool(self, width=0.5):

		if self.redraw:
			glNewList(self.listids, GL_COMPILE)
			self.drawBand(width)
			self.draw()
		#	self.drawTube(0.5,16)
			glEndList()
			glCallList(self.listids)
			self.redraw = False
		else:
			glCallList(self.listids)
